Muffin Tin Orbitals of Arbitrary Order 

O.K. Andersen and T. Saha-Dasgupta 
Max- Planck- Institut fur Festkorperforschung, Postfach 800665, D- 70506 Stuttgart, Germany 

(February 1, 2008) 

Abstract 

We have derived orbital basis sets from scattering theory. They are ex- 
pressed as polynomial approximations to the energy dependence of a set of 
partial waves, in quantized form. The corresponding matrices, as well as the 
Hamiltonian and overlap matrices, are specified by the values on the energy 
mesh of the screened resolvent and its first energy derivative. These orbitals 
are a generalization of the 3rd-generation linear MTOs and should be useful 
for electronic-structure calculations in general. 
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For electrons in condensed matter, it is often desirable to express the one-electron wave 
functions, (r) , with energies, in a certain range in terms of a minimal set of energy- 
independent orbitals, xrl (r) . Here, R labels sites and L the local symmetry (e.g. L=lm). 

The simplest example of such an orbital is the Wannier function, \ (r — R) , for an 
isolated band. A more realistic example is illustrated in Fig. |I], the conduction-band orbital 
of a cuprate high-temperature superconductor. This orbital is centered on Cu, has anti- 
bonding O x p x -Cud x 2_ y 2 -OyPy character, and extends beyond the 3rd-nearest neighbor 
atoms. Its Bloch sum describes a tight-binding (TB) band: ~ (s) — It (cos k x + cos k y ) + 
At' cos k x cos k y — It" (cos 2k x + cos 2k y ) . This orbital is the starting point for descriptions of 
the low-energy physics of the cuprates. Its is not a Wannier function. First of all because 
the conduction band is merely one partner of a bonding, non-bonding, anti-bonding triple 
with nearly degenerate Cu d and O p levels so that the three bands nearly stick together at 
Ep^Ed with a cone-like behavior at the centre of the zone. As a result, the true Wannier 
function of the anti-bonding band has very long range, but since E p ~Ed is 2-3 eV below the 
Fermi level, the low-energy physics is hardly influenced by this. The second reason why the 
orbital of interest cannot be a Wannier function, is that the conduction band is crossed by, 
or has avoided crossings with other bands (Fig.|2|). Since this occurs an eV below Ef, this, 
too, is irrelevant for the low-energy physics, which should therefore be described using an 
orbital which yields correct wave functions at and near Ep and has errors oc (e$ — Ep) N+1 . 
The wider the energy range described correctly by this orbital, i.e. the higher the N, the 
longer its spatial range. 

We have found a general method, the NMTO method, by which for instance this kind 
of orbital can be obtained [|]]. What Fig. [I] shows is in fact a muffin-tin orbital (MTO) with 
N—l, obtained from a density-functional (DF-LDA) NMTO calculation. This method has 
recently enabled us to compute how the hopping integrals t, t', and t" are influenced by 
chemical and structural factors, and it has proved successful for computing t\\ and t± for the 
ladder cuprates without resort to the common, but dubious procedure of fitting to guessed 
TB bands [§. 
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In Fig. |2| we demonstrate that a single MTO of sufficiently high N is capable of describing 
the entire conduction band, including its cone-like feature as well as smooth interpolations 
across avoided crossings: The dotted band was obtained variationally using an MTO with 
N=3, thus yielding band-errors of order 2(iV + 1)=8. This figure also demonstrates that 
one may use a discrete mesh of energies, eo, ejv, to construct the MTO, which then has 
errors oc (e, — eo) ... (e* — ejv) • This is analogous to using Lagrange or Newton interpolation 
instead of Taylor expansion, and is far more practical. The band obtained variationally has 
errors oc (e* - e ) 2 ... (e* - e N f . 

For some purposes, it is better to use a larger set of more localized orbitals. For instance, 
in order to understand the microscopic origins of t, t', and t", we used a set with Cu d x 2_ y 2, 
Op x , Op y , and Cus, obtained by upfolding through a screening transformation |L],|]||. 

Materials with many bands and strong correlations are being studied intensively. The 
first step of a quantitative description is a one-electron mean-field theory requiring a basis, 
flexible enough to give individual orbitals desired properties. For this, NMTOs are uniquely 
suited. 

As an example of a minimal set spanning all states in a wide energy range, let us consider 
the LDA valence and conduction bands for GaAs, 18 of which fall in the range between - 
15 and +20 eV. With a G&sp 3 d 5 As sp 3 d 5 f 7 basis of merely 25 N=2 MTOs per GaAs, and 
mesh points at -15, 0, and 10 eV, we obtained a variational band structure, which only 
above +15 eV yielded errors as large as 0.1 eV. Even for this 35eV-range, which includes 
the Ga 3c? semi-core band at -15 eV, no principal quantum numbers were needed. To most 
practitioners, this is surprising result. NMTOs should be useful for computing excited-state 
properties with the GW method 0. 

For ground-state properties, only the Ga 3<i and the valence bands must be described. 
Using the minimal Ga sp 3 d 5 As sp 3 MTO set, we find accuracies in the sum of the one- 
electron energies of 50 and 5 meV per GaAs for respectively N—l and N=2 |T|]. This is 
highly satisfactory and opens the way for accurate and efficient DF-calculations, for instance 
for large systems using techniques where the computation increases merely linearly with the 
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size of the system. Hitherto, this has only been possible with less accurate or geometry- 
restricted methods such as semi-empirical TB, screened LMTO-ASA [|J, or screened 
multiple-scattering theory |Q. 

The LMTOs of the 1st- and 2nd-generations were expressed in terms of partial waves, 
ifui (e , r R ) Y L {f R ) , and their energy derivatives, (pjy (e , r R ) Y L (r R ) , truncated outside the 
atomic spheres (r R = |r — R|). Everything else was neglected in the atomic-spheres approx- 
imation (ASA), which then gave rise to a simple formalism and fast computation. The 
3rd-generation || succeeds in making this formalism valid for overlapping MT potentials, 
^ ( r ) = 12r v R { r n) 5 to first order in the overlap of the v's, thus making the ASA superfluous. 
This is accomplished by attaching tails of screened spherical waves with the proper energy to 
the partial waves. The resulting set of kinked partial waves, evaluated on the energy mesh, 
is what the NMTO set is expressed in terms of: 

This may be considered as a polynomial approximation to the energy dependence of the 
partial- wave set, in quantized form. In the following, we derive the expressions for the 
Lagrange matrices, Ln , and the NMTO Hamiltonian and overlap matrices, starting out 
from the conceptually simplest way of solving Schrodinger's equation, namely by matching 
of partial solutions. Our formalism should prove useful also in other contexts. 

We consider the case where the wave functions (r) are solutions of a Schrodinger 
equation with a MT potential, Ti^i (r) = [—A + V (r)] ^ (r) = e^i (r) . For simplicity, we 
first assume that the MT wells do not overlap and have ranges, a R . At the end, definitions 
will be modified in such a way that the formalism holds also for overlapping wells. The a's 
will be hard-sphere radii which define the screening and, hence, the shape of the orbitals. 

Kinked partial waves 0. -Inside a MT sphere, the partial solutions factorize into 
energy-dependent radial functions, <pm (e, r R ) , and angular functions. In the interstitial, 
we use screened spherical waves, which are defined as those solutions of the wave equation, 
(A + e) ip R L (e, r) = 0, which satisfy the homogeneous boundary condition that the projec- 
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tion ofipRL (s, r) onto S (tri — or/) Yu {jrj) be 5rri5ll>- In fact, only those solutions with RL 
corresponding to the so-called active channels will be used (in Fig. [I], the central Cu d x 2_ y 2), 
and only the projections onto other active channels will vanish (all non-central Cu d x 2_ y 2 
projections). The projection of i/)rl i £ , r ) onto an inactive channel (all other than d x 2_ y 2 
on any Cu-sphere) satisfies the boundary condition that its radial logarithmic derivative 
equals that of the radial Schrb ding er-solution. The kinked partial wave, <J)rl (e, r) , is now 
fm { £ , t r) Yl (j'r) inside its own sphere and for its own angular momentum, it is ipRL (e, r) 
in the interstitial region, and inside the sphere at R', it vanishes for any other (R'L'^RL) 
active channel, but is proportional to (pR>i> (e,rR>) Yjj {tri) for an inactive channel. As a re- 
sult, with the normalization ipRi (e, a#) =1, the kinked partial wave is a continuous solution 
of Schrodinger's equation with energy e. But since it has kinks at the spheres in the active 
channels, it is not a wave function. 

The solid curve in the left-hand part of Fig.|3| shows the Si p x=y=z kinked partial wave for 
e in the middle of the valence band and for r along the [lll]-line in the diamond structure 
from the central Si atom, through the nearest Si neighbor, and half-way into the back-bond 
void. The other curves will be explained when we come to consider potential overlap. The 
kinks at the a-spheres (chosen smaller than touching) are clearly seen. Since this kinked 
partial wave is designed for use in a minimal sp 3 -basis, only the Si s and p waves were 
chosen as active. The inactive waves must therefore be provided by the tails of the kinked 
partial waves centered at the neighbors, and this is the reason for the strong Si <i-character 
seen inside the nearest-neighbor sphere. Had we been willing to keep Si <i-orbitals in the 
basis, the Si <i-channels would have been active so that only waves with l>2 would have 
remained inside the neighbor spheres, whereby the kinked partial wave would have been 
more localized. Hence, the price for a smaller kinked-partial wave basis, is longer spatial 
range and a stronger energy dependence. 

The element Kr>l',rl (e) of the Hermitian kink matrix is defined as the kink of 4>rl (e, r) 
at the a#/-sphere, projected onto Yl> (rn>) I^ri- Hence, it specifies how the Hamiltonian 
operates on the set of kinked partial waves: 
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{H - e) <\> RL (e, r) = [-A + V (r) - e] <\> RL (e, r) = (2) 
- ^ 5 (rjy - or/) F L / (r R *) K R , L ^ RL (e) . 

Although an individual kinked partial wave is not a wave function, any smooth linear combi- 
nation, J2 R l (e, r) crl,% > is. Schrodinger's equation may therefore be formulated as the 
matching- or kink-cancellation condition: J2 R l^-R'L',rl (Si) c R L,i = for all R'L', which is 
a set of homogeneous linear equations, equivalent with the KKR equations §J. Here, the 
indices run only over active channels. Since the kink-matrix is expensive to compute, it is 
not efficient to find a one-electron energy from: det \K (£j)| =0, and then solve the linear 
equations for the corresponding c R L,i- Rather, we construct a basis set, x ( r ) > with the 
property that it spans any wave function, (r) , with an energy in the neighborhood of 
iV+1 chosen energies, €o,...,€n, to within an error oc (e» — 6q) ... (e* — ejv), and then solve 
the generalized eigenvalue problem, 

^2 RL (x i £l\'H-^\x { S)baL,i = for all R'L', (3) 

resulting from the Raleigh-Ritz variational principle. 

MTOs. -Since all wave functions with Si=e may be expressed as: YIfll ^RL ( £ i r ) c RL,i, 
the MTOs with iV=0 are simply the kinked partial waves at the chosen en- 
ergy: Xrl ( r ) —^RL (^q, r) . The Hamiltonian and overlap matrices are respectively 
<X (0) |^- e o|x (0) > = -K(e ) and <x (0) | X (0) > = K (e ) , as may be found from Eq. (|) 
and the normalization chosen. Here, ~=d/de. In order to find the MTOs with iV>0, we 
first define a Green matrix: G (e) = K (e)" 1 , and then, by an equation of the usual type: 
(H — e) 7i?L (e, r) = —5 (r R — a R ) Yl (r R ), a Green function, j R l (e, r) , which has one of its 
spatial variables confined to the a-spheres, i.e. r'— >RL. Considered a function of r, this 
confined Green function is a solution with energy e of the Schrodinger equation, except at 
its own sphere and for its own angular momentum, where it has a kink of size unity. This 
kink becomes negligible when e is close to a one-electron energy, because the Green function 
has a pole there. Eq. (g) shows that 7 (e, r) = (e, r) G (e) . (Here and in the following, 
lower-case letters, such as 7 and 0, denote vectors, and upper-case letters, such as K and 

6 



G, denote matrices; e, e, RL, and N are numbers, though). The confined Green function is 
thus factorized into a Green matrix G (e) which has the full energy dependence, and a vector 
of functions <fi (e, r) which has the full spatial dependence and a weak energy dependence. 
(The energy windows we consider are limited in size by the requirement that <Prl (e, r) and 
4>rl (s'j r ) cannot be orthogonal). Finally, we want to factorize the r and e-dependences 
completely and, hence, to approximate the confined Green function by x ( r ) G (e) : We 
note that subtracting from the Green function a function which is analytical in energy, 
(e, r) G (e) —oj^ (e, r) = x { £ , r ) G ( £ ) > produces an equally good Green function in the 
sense that both yield the same Schrodinger-equation solutions. If we can therefore determine 
the vector of analytical functions, u>^ (e, r) , in such a way that each Xr2 ( £ , r ) takes the 
same value, Xr2 ( r ) j a ^ a ^ mes h points, then Xrl ( £ > r ) = Xrl ( r ) + O ((e — e ) .. (e — ejy)) • 
Hence, x (A ° (r) is the set of NMTOs. Now, since x (N) (e , r) = ... =x (A ° (e^, r), the iVth di- 
vided difference of (s, r ) G (s) equals ( r ) times the Nth. divided difference of G (e) . 
Moreover, if we let (e, r) be a polynomial in energy of (iV-l)st degree, its iVth divided 
difference on the mesh, A N u^ (r) /A [0...7V] , will vanish. We have therefore found the 
following solution: 



X W (r) 



A N (/> (r) G 
A[0...iV] 



A N G 
A]oZ/V] 



(4) 



[JV-1,JV] 

■■♦lio^-*)"^-^- < 5 > 

for the NMTO set. Since the kinks, (H-e)<t> (e, r) G (e) , are independent of e, NMTOs 
with iV>0 are smooth. By use of the well-known expression for a divided difference: 

A"*(r)G ^ «!(£„, r)G(e„) 



1 n=0 1 lm=0,^n V c 



n=0 1 lm=0,^n V L n e m 



we finally obtain the expressions for the Lagrange matrices in Eq. (0) and the energy matrices 
in Eq. <g): E {M ^= (A M eG/A [0..M]) (A A/ G/A [0..M])' 1 , in terms of the values of the Green 
matrix on the energy mesh. 



The NMTO set may thus be thought of as a 'quantized' Lagrange interpolation of the 
kinked partial- wave set, where the weights are matrices rather than iVth-degree scalar poly- 
nomials in energy. Similarly, Eq. (^) may be interpreted as a 'quantized' Newton interpo- 
lation with the energies substituted by matrices. If the mesh is condensed, Newton in- 
terpolation becomes Taylor expansion: A </>/A [0...iV] — > (1/iV!) d N <f)/de N . The form ([S]) 
expresses the NMTO as a kinked partial wave at the same site and with the same angular 
momentum, plus a smoothing cloud of energy-derivative functions centered at all sites and 
with all angular momenta. In the right-hand part of Fig. [3], the solid curve is the MTO 
with N—l, and the dashed curve is the MTO with iV=0 shown also in the left-hand part. 
Here again, longer spatial range is the price for spanning the wave functions in a wider 
energy range. The increase of range and smoothness with N follows from the relation: 
(Tt — ejv) x ( r ) — X ( r ) — Cat) , which also shows that the £"s are transfer ma- 

trices between MTO sets of different order. Linear transformations of the kinked partial 
waves, 4> { e i r ) =4> ( £ ) r ) T (e) , change the NMTOs, but not the Hilbert space spanned by 
them |[J. This may be used to generate nearly orthonormal representations where the E's 
are Hamiltonians and where (x 1 I X ) = 1 f° r 1 < M < iV. 

The expressions for the Hamiltonian and overlap matrices needed in (^) may be worked 
out and given as 0: 

X (JV) k-W|x (JV) >A4n= ( 6 ) 



A[0...iV] XA 1 IA ' A[0...N] 

A 2NQ A 2N+1 G 



(e - e Nj 



A[[0..N-1]N] v ' A[[0...iV]]' 
A M+N+1 G/A [[0..M] N] is the (M+iV+l)st derivative of that polynomial of degree M+N+l 
which takes the values G (eo) , G (ejv) at the iV+1 mesh points and, at the first M+l 
points, also the values G (eo) , .., G (6m) of the energy-derivatives. The one-electron energies 
are 'ratios' of energy derivatives of such 'Hermite interpolations' of G (e) , which itself has 
poles inside the mesh. 

Having seen that the formalism is expressed 

in terms of one matrix, e.g. K (e) = \ £ — ^1 X^) = G , let us indicate how this 
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is generated P,|iO|: The elements of the bare KKR structure matrix ||, B R , L , RL (e) 



Y lP ±™~ l¥l '~ l "C LV i»ivrn» («|R-R'|)y^, (r^R') for R^R', and =0 for R=R', spec- 
ify how the spherical waves, n\ (ktr) Yl (tr) , are expanded in regular spherical waves, 
jv (^r')Yl' (Pr') ■ The corresponding expansions of the screened spherical waves are now 
specified by a screened structure matrix, defined via: B a (ey 1 = B° (sy 1 + K~ x tana (e) , 
and obtained by matrix inversion of B° (e) +k cot a (e) . Here, k cot a (e) is a diagonal matrix 
with am (s) being the hard-sphere phase shift, tan a ri (e) =ji (kor) jn\ (kclr) , if the chan- 
nel is active, and the true phase shift, T]ri (e) , if the channel is inactive. B a (e) has short 
spatial range for energies well below the 'hard-sphere continuum,' as defined by the division 
into active and inactive channels and the choice of a-radii for the former. The kink matrix 
is finally: K (e) = — [nn {ko)\~ [B a (e) + K cot 77° (e)] [nn (na)]^ 1 , where r] a (e) is the phase 
shift in the medium of hard a-spheres: tan i] a (e) = tan rj (e) — tan a (e) . B a (e) contains the 
essence of the hopping integrals, whose dependence on the local environment enters through 
the screening. 

When the potentials overlap, we need to redefine the kinked partial waves as illustrated in 
Fig.0: <f) RL (e, r) = [y? ffl (e, r R ) -(p° m (e, r R )} Y L (r R ) +i) RL (e, r) . Here, ip (e, r) (dot-dashed) is 
the radial solution for the central MT-well, which now extends to s (> a). (f° (e,r) (dotted) 
is the phase-shifted wave proceeding smoothly inwards from s to the central a-sphere, where 
it is matched with a kink to the screened spherical wave ip (dashed). It is easily shown that, 



with this modification, the formalism holds to first order in the potential-overlap [|I|,|5|,|10 
In practice, this means that radial overlaps of up to 30% may be treated without changes, 
and that overlaps as large as in Fig.|3|, may be treated by adding a simple kinetic-energy 
correction [|T|j5| JT0|JTTJ| . This should make the use of empty spheres superflous and open the 
way for efficient DF-molecular-dynamics calculations. The a-radii now specify the screening, 
with a default value which is 80% of the atomic or ionic radius, and for semi-core states, the 
core radius. 

In conclusion, we have solved the long-standing problem of deriving useful, minimal sets 
of short-ranged orbitals from scattering theory. Into a calculation enters: (1) The phase 



shifts of the potential wells. (2) A choice of which orbitals to include in the set, the so- 
called active channels. (3) For these, a choice of screening radii, or^, to control the orbital 
ranges. (4) An energy mesh on which the set will provide exact solutions. These MTOs 
have significant advantages over those used in the past. 
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FIGURES 

FIG. 1. The Cu d^^-like LMTO, which describes the (LDA) conduction band of HgE^CuOzt, 
plotted in the Cu02 plane. Cu and O sites are marked by respectively + and 

FIG. 2. Band structure of CaCu02 with a 7°-buckle, calculated in the LDA with a single Bloch 
Cu d x 2_ y 2 CMTO (dotted) compared with the full band structure (solid). 

FIG. 3. Si p x=y=z kinked partial wave (KPW), its constituents ip,ip°, and ip, and the LMTO. 
No empty spheres were used, s is the range of the central potential well. 



12 



